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^ INTRODUCTION 



ABSTRACT 

Weak gravitational lensing on a cosmological scales can provide strong constraints both on the nature 
of dark matter and the dark energy equation of state. Most current weak lensing studies are restricted 
to (two-dimensional) projections, but tomographic studies with photometric redshifts have started, and 
future surveys offer the possibility of probing the evolution of structure with redshift. In future we will be 
able to probe the growth of structure in 3D and put tighter constraints on cosmological models than can 
be achieved by the use of galaxy redshift surveys alone. Earlier studies in this direction focused mainly on 
evolution of the 3D power spectrum, but extension to higher-order statistics can lift degeneracies as well 
as providing information on primordial non-gaussianity. We present analytical results for specific higher- 
order descriptors, the bispectrum and trispectrum, as well as collapsed multi-point statistics derived from 
them, i.e. cumulant correlators. We also compute quantities we call the power spectra associated with the 
bispectrum and trispectrum, the Fourier transforms of the well-known cumulant correlators. We compute 
the redshift dependence of these objects and study their performance in the presence of realistic noise 
and photometric redshift errors. 

Key words: : Cosmology- Weak Lensing Surveys- Large-Scale Structure of Universe - Methods: ana- 
lytical, statistical, numerical 



tUVnil very recently, the best information about the power spectrum of cosmological density perturbations has been obtained from large-scale galaxy 
-surveys and Cosmic Microwave Background (CMB) observations. However, galaxy surveys only probe directly the clustering of luminous matter, while 
>^MB observations mainly explore the power spectrum at a very early linear stage of their evolution. Weak gravitational lensing stud ies provide a 
J6»mplementary approach for probing the cosmological power spectrum at modest redshift in an unbiased way; for a recent review, see Munshi et al. 
( 2008). Weak lensing is a relatively young subject; the first measuremen ts were published within the last decade dBeacon. Refregier & Ellis f 20001 ; 
IWittman et alll200ol ; lKaiser. Wilson & Luppinoll2000l ;[Waerbeke et al 2000). Since then rapid progress has been made on analytical modelling, technical 
specification and the control of systematics. Over the course of the next few years, photometric redshift surveys are going to be increasingly prevalent. 
Such deep imaging surveys combined with resulting photometric redshift information will mean that there will be a considerable scope for using weak 
lensing studies to map the dark matter in the universe in three dimensions. 

Ongoing and future weak lensing surveys such as the CFHT legacy survejQ, Pan-STARRS0, the Dark Energy Survey, and further in the future, the 
Large Synoptic Survey Telescope 0, JDEM and Euclid will provide a wealth of information in terms of mapping the distribution of mass and energy in 
the universe. With these surveys, like the CMB, the study of weak lensing is entering a golden age. However to achieve its full scientific potential, control 
is needed of systematics such as arise from shape measurement errors, photometric redshift errors, an d intrinsic alignm ents. 

Th e use of photometric redshifts to study weak lensing in three dimensions was introduced by iHeavensI |2003). It was later developed by many 
authors dHeavens et alj|200ol ; iHeavens et all 2006 ; iHeavens, Kitching & Verde]|2007l; ICastro et alll2005h , and was shown to be a vital tool in constraining 
dark energy equation of state dHeavens et all200q) . neutrino mass dKitching et alj2008h and many other possibilities. While the traditional approach deals 
with projected surveys or with tomographic information, there has been substa ntial recent progress in the d evelopment of studying weak lensing in 3D. 

Early analytical work in weak lensing mainly adopted a 2D approach djain. Sel iak & White 2000), due to the lack of an y photometric redshift 
information about the source galaxi es. It is also interesting to note that most of these studies also employed a flat-sky approach dMunshi & Jainll200l[ 



information about the source galaxies. It is also interesting to note that most of these studies also employed a flat-sky approach (Munshi & Jain 2001; 
lMunshill2000l ; lMunshi & J ain 2000), as the first generations of weak lensing surveys mainly focused on small patches of the sky dMunshi et al.ll2008h . 



1 http://www.cflit.hawaii.edu/Sciences/CFHTLS/ 

2 http://pan-staiTS.ifa.hawaii.edu/ 

A http://www.lsst.org/lsst_home.shtml 
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T hese works made analytical predictions for lower-order moments as well as the entire probability distribution function of convergence field k or shear 
7 dMunshi & Jainll200ll ; IValageasl200ol : lMunshi & Valageasll2005l ; IValageas. Barber. & Munshill2004l ; IValageas. Munshi & Barberjl2005l). A tomographic 
("2.5D")approach has also been developed , wherein the sources are divided into a few redshift slices dTakada & W hite 20031: iTakada & Jain] 120041 : 
iMassev et alll2007l : ISchrabback et al.| [2009) and these slices are then analyzed jointly, essentiall y using the 2D approach but keeping the information 
regarding the correlation among these redshift slices. A notable exception however in this trend was Stebbins(1996) who developed the analysis techniques 
for weak lensing surveys covering the entire sky. In recent years there has been a lot of interest in developing analysis tools and predicting the cosmological 
impact of future generations of weak lensing surveys with large sky coverage which are naturally analyzed in the spherica l harmonic dom ain. I n this paper, 
we pre sent a very general study of 3D Weak lensing beyond the power spectrum. Extending the formalism developed in [Heaven] {2003) and lCastro et all 
(2005) we use higher-order statistics to probe non-Gaussianity present in primordial anisotropy as well as that induced by gravity. 

While power spectrum analysis does provide the bulk of the information regarding background cosmology, higher-ord e r statistics are useful t o lift 
degen e racies, allowing determination of Q m and as independently - see, e.g.. IBemardeau. V an Waerbeke & Mellierl dl997h : |jain & Seliakl ( ll997l) : lHul 
( 1999); Schneider et alt d2002l) ;|Takada & Jain] d2003h . Some of these studies also carried out tomographic analysis using the bispectrum. Detailed Fisher 
matrix analysis found that the level of accuracy with which various cosmological parameters including the dark energy equation of state parameters 
can be enhanced considerably by using bispectrum data in combination with power spectrum measurements. Higher-order studies also are important fo r 
evaluating scatter in lower-order estimates; e.g. the trispectrum is important for computing error bars in power spectrum estimates dTakada & Jaiil2 009). 
Most studies involving higher-order correlations however mainly concentrate on projected convergence. The aim of this paper is to extend higher-order 
analysis to 3D by taking into account the radial distance as inferred through photometric redshift information. 

Modelling of the underlying mass distribution is necessary for predictions of weak lensing multispectra . In earlier studies, the hierarchical ansatz 
was found to be very useful in modelling higher-order statistics of weak lensing ob s ervables dFrvlll984l; ISchaefferjll984l; lBernardeau~ & Schaeffer 
1992l:ISzapudi & Szalavll993|.ll997l:lMunshi et alll999l : lMunshi. Coles & Melottll999al lb1: lMunshi. Melott & Colesll999l ; IColes. Melott & Munshill 19991 : 



Munshi & C oles 2000, 2002, 2003). The hierarchical ansatz models the higher-order correlations as a hierarchy, with higher-order correlation functions 



being expressed as products of correlation functions. The diagrammatic representation of these expressions resembles perturbative models of the corre- 
lation hierarchy which typically develops with the onset of gravitational clustering in collisionless media. The amplitudes of these "Feynman diagrams" 
are of course different in the perturbative regime and in the quasilinear regime. Various hierarchical ansatze differ in the way they assign amplitudes to 
diagrams with different topologies. We will employ the most generic hierarchical ansatz in modelling the underlying mass distribution, and the method 
can be modified to take into account any other specific forms of correlation hierarchy in a rel ative l y simp l e way. 

Higher-order correlation functions have also been detected observationally dHuil 1 19991 : IBemardeau. Van Waerbeke & Mellier] 1 19971 ; 
IBemardeau. Mellier & Van W aerbeke 2002). As expected though these studies are more difficult than two-point as they can be domi nated by noise. There 
have been several studies in this di rection which focuses mainly on projected surveys as well as using tomographic information (Hu 1999; T akada & Jainl 
l2004l2003l ; ISemboloni et alll2008h . Typically one-point cumulants or lower-order moments are employed to compress information in higher-order cor- 
relation functions into a single number. This is indeed due to the related gain in signal-to-noise. A full analysis of multi-point correlation function (or 
their respective Fourier transform s, the multispectra) is relatively difficult because of the low signal-to-noise ratio of individual modes. In their recent 
study iMunshi & Heavens! {2009) used an intermediate option: they found that a better approach, one that is even optimal in certain cases, is to use the 
power spectra associated with the various multi-spectra. These objects combine various individual modes of the multi-spectra in a specific way and can 
be computed from numerical simulations or observed data relatively easily. While their work was motivated by cosmological studies of non-Gaussianity 
involving the CMB, here we generalize it to 3D weak lensing studies. We primarily focus here on convergence studies but the results can be generalized 
to shear statistics. We focus on three- and four-point statistics, but the formalism is general. We develop the analysis tools and provide results both 
in all-sky limit using a harmonics treatment (useful for future surveys with large sky coverage) as well as using a patch-sky approach using flat-sky 
Fourier transforms (for surveys with relatively small sky coverage). In addition to full analytical results we also provide results using the extended Limber 
approximation which can drastically reduce the computational cost with very little loss of accuracy at high wavenumbers. 

The paper is arranged as follows: In §2 we discuss the basic formalism of 3D weak lensing and how it can be used to estimate the power spectrum 
in 3D. It introduces the notations which will be used in the following sections. In §3 we introduce the models describing higher-order clustering and their 
associated hierarchical amplitudes which are then used to construct a model for the bispectrum and higher-order multispectra in the nonlinear regime. In 
§4 we focus on representation of bispectrum and trispectrum respectively in various coordinate systems. In §5 we relate observed convergence statistics 
with the underlying statistics of mass distribution. Sections §6 and §7 are devoted to development of the skew spectrum (3-point) and the kurt-spectrum 
(4-point). In §7 we develop the for malism for surveys with small sky coverage and §8 is reserved for discussion of results and future prospects. Throughout 
we will borrow the notations from lCastro et all d2005l) wherever possible. We have analysed the effect of photometric redshift errors as an appendix. 

The general formalism developed in this paper will have applicability in other areas of cosmology, where 3D information can be used effectively, 
including future generations of 21cm surveys as well as near all-sky redshift surveys. 



2 WEAK LENSING IN 3D: ANALYSIS OF THE 3D POWER SPECTRUM 

In this section we introduce the formalism of weak le nsing in 3D. It w ill also serve to introduce the notations which we will u se in the later se ctions. 
The formalism for 3D w eak lensing was de veloped by [Heavens] d2003l) and was developed further by several authors including ICastro et al d2005h . We 
will review results from ICastro et all {2005) for power spectrum analysis here which will be useful for introducing t he notations and for t he fur ther 
developments presented in later sections. We will also use the extended Limber approximation (to be introduced later; iLo Verde & Afshordl d2008l) ) to 
simplify results obtained in Castro et al (20051). 

In this work we will mainly be concerned with the convergence field. Analytical extensions to deal with shear fields will be dealt with elsewhere. 
We introduce $(r, 0, ip) as the 3D gravitational potential at a 3D position r, 9, <p, and 0(r) the lensing potential. The radial distance r(t) is related to the 
Hubble expansion parameter H(t) = a/a by r(z) — c J dz' /H(z'). The Hubble parameter is sensitive to the contents of the Universe thereby making 
weak lensing a useful probe to study dark energy. The line of sight integral relating the two potentials can be written as dKaiserlll992l) : 
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Deriv ation of the above expression assumes the Born approximation ( iBernardeau. Van Waerbeke & Mellienl 19971: ISchneider et aj|2002l : IWaerbeke et al 
2002), which evaluates the line of sight integral along the unperturbed photon trajectory. Note that the lensing potential cj>(r, fi) is radially dependent. 
Throughout these papers we denote vectors in bold letters, c denotes the speed of light, r — r(t) is the comoving distance to the source at a given instance 
of time t from the observer who is situated at the origin (r = 0). Depending on the background cosmology /ic(r) can be sin r, r, sinh r for a closed 
(K = 1), flat (K = 0) or open (K = — 1) universes respectively. Our convention for the Fourier transform for the 3D fields closely resemble that of 
ICastroetallj2005h . The eigenfunctions of the Laplacian operator in flat space when expressed in spherical coordinates turn out to be a product of spherical 
Bessel functions ji(kr) in the radial direction and the spherical harmonics on the surface of a unit sphere i.e. Y; m (fi) = Yi m (6, (f>). The eigenfunctions 
ji(kr)Yi m (8, cj>) are associated with eigenvalues — k 2 . The eigendecomposition and its inverse transformation can be expressed as: 

j ~ p l~~ p oo m — l 

$lm(k) = J^ / d^4(r)*ji(to0*im(ft)> $(r) = w|; / kdk^ *h»(*)jl(*r)y* m (ft). (2) 

" J " * l=Q m =-l 

The choice for this eigendecomposition is determined by various factors. First it can deal with large areal sky coverage, and secondly as the lensing is 
related to gravitational potentials, the expansion allows us easily to expre ss the coefficients of expansion of the convergence (or shear) in terms of the 
expansion of the density field through the Poisson equation Heav ens! d2003l) . &im(k) here is the spherical harmonic decomposition of $(r), and similarly 
for <f>(r). The orthogonality properties for the harmonic modes for the 3D potential can be used to define the 3D all sky power spectra by the following 
expressions. 

(*h,(*)*iw(*')) = C?*(k)5 1D (k + k')8%8K m ,; (<f>i m (k)4>i>m>(k')) = C^(k)8 1D (k + fe')4<*w • (3) 

Here the power spectrum C** represents the 3D power spectrum associated with the 3D potential field $; m (fc) and 8„d and 8 K represent the n- 
dimensional Dirac and Kronecker delta functions respectively. The corresponding all-sky power spectra for the lensing potential <f> is denoted by Cf^. In 
comoving coordinates we can write: 

A<f(r) = ^-n m Hl5{r); $ im (fc; r) = -^-l^Q^S^k; r) = -^S lm (k; r). (4) 

Here a(z) = 1/(1 + z) is the scale factor at redshift z, f2 m is the total matter density at z = 0, and Ho is the Hubble constant today. 8i m {k; r) is 
the eigendecompo sition of S(r). Wh en appearing after the semi-colon, the r dependence (e.g. of $; m (fc;r)) really expresses the time-dependence of 
the potentials; see lCastro et all ( 120051) for a discussion of the subtleties of this. Any model we assume for describing non-linear growth of perturbations 
for S(r) will thus have direct impact on statistics of observed weak lensing convergence « as it depends on <5(r) through i ts dependence on $(r). The 
harmonic decomposition of the lensing potential and the 3D gravitational potential are related by the following expression (Castro et al 20051): 



»W — — o / dk'k' I rdrji(kr) j dr' 



ji(k'r')$ lm (k';r'). (5) 

We ignore complications of varying radial selection function for the sources, and distance errors in the main text for clarity. These are considered in 
Appendix A. 

Using the relation which relates the density and the gravitational potential as well as expressing the convergence field in terms of the lensing potential 
Kim(fc) = — + l)<f>i m (k)', we can express the convergence coefficients in terms of those describing the density field. This is important because we 
can then relate the statistics of the density field with that of the convergence field directly which is potentially observable: 



. /, / /\ 8im(k ; v ) 



We have absorbed the cosmological constants in the constant A = — 3£Im Hq /2 We will also introduce the quantity Ii(ki, k), which will be useful in 
displaying the future results: 

h{ki,k) = ki( drrji{ki,r) [ dr' ( T —^- J j,(kr')y/P**(k; r>) (7) 



We can now write down the power spectrum associated with the lensing potential in a more compact form. Next we relate the convergence power 
spectra in terms of the Cf using their relationship in the harmonic domain. 

1fi4 2 f°° 1 
C?*(fci,fca) = -=-=- / fc 3 I,(fci,fc)I,(*a,fc)dfci Cr(*i,fca) = 7 I a (i + l) a Cf*(fci,fca). (8) 
*" c Jo 4 

In deriving the above expression it was assumed that gravitational potential power spectrum can be accurately approximated by P ~ 
/ P**(fc; r)P**(fc; r 1 ). For a detailed description and range of validity see ICastro et a il d2005l) . Clearly the analysis outlined above follows three 
different steps. First we relate the lensing potential 4>i m (k; r), or equivalently Ki m (k; r), to the 3D potential $; m (fc; r). Next the statistics of $; m (fc; r) 
are used to make concrete predictions about the st atistics of K; m (fc; r ). However an intermediate step is required to connect the 3D Fourier decomposition 
<5(k) and its harmonic counterpart 8i m (k; r). See jCastro etalll2005l) for details. 

Throughout this paper, the analysis will rely on various assumptions. For an arbitrary field Vl/(r; r) which is assumed isotropic and homogeneous with 
spherical harmonics decomposition ^i m (k) can be characterized by a power spectrum Ci(k; r). It is important to realise that r, the comoving distance, 
also plays the dual role of cosmic epoch, and we retain r in harmonic representations to label the cosmic epoch. The cross-power spectrum related 
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to an arbitrary 3D field at tw o different radial d istances (redshifts) Ci(k;r,r') will be expressed as {ty i m (k; r)^*,^ (k' ; r')) — G'i(k;r,r')SiD(k — 
k')5^i8^ m i . It was shown bv lCastro et all 120051) that Ci (k, r) is simply the 3D power spectrum P(k; r), Ci (k; r, r') — P(k; r, r'). For the derivation 
we need to expand the Fourier decomposition of ^(r; r) and exploit the fact that (^(k; r)^'(k; r')} = (27r) 3 P(&; r)(53L)(k — k'). In our present analysis 
we will focus on extending these results to higher-order statistics. We model the underlying statistics of the density field by using non-perturbative results 
and relate these to the statistics of projected field such as convergence. We also introduce power spectra related to multispectra to effectively compress 
the information content. The results are presented both in harmonic space as well as in Fourier domain using the Fourier approximation. 



2.1 3D Convergence Power Spectrum Using the Limber Approximation 

Computations of higher-order multi-spectra are often difficult, given the mul ti-dimensional integrals involved, which often make numerical computations 
time consuming if not prohibitive. The Limber approximation dLimbed 19541) or its generalization to Fourier space is often used to simplify the evaluation 
numerically by reducing the dimensionality of the integrals. Typically implementation of the Limber approximation is valid at small angular separations 
which correspond to large multipole moments I in harmonic domain. It also requires smooth variations of the integrand compared t o the Bessel functions 
of rele vant I. For a detailed description of various issues and calculations of next order correction terms see a recent discussion by Lo Verd e & Afshordil 

fcoosh . 

We start by the following expression Eq.(|5} and Eq.([8} from the previous section: 



16A 2 



2„2 



7T Z C- 



k dkli(ki, k)Ii(k2, k) 



ki / dr a r a ji(kir a ) / dr' a 



k 2 ] dnr b ji(kir h ) / dr' b 



/ dr' b 


n - r' b 


lo 





x / fc 2 dfc v /p**(fc;^)P**(fc;<)j i (fcr a )j,(fcr 6 ) 



(9) 

We now use the Limber approximation to simplify the k integral which produces a 8iD(r' a — r' b ) function. Integrating out r' b with the help of the 
delta function and renaming the dummy variable r' a to r we can finally write: 



cr(ki,k 2 ) 



7vl 2 (l + lfA 2 



Ii(ri,r 2 ) = 



4 
16 



fei k 2 



dr 1 j l (kir 1 ) 



dr' 



r 2 



dr 2 ji(k 2 r 2 )Ii(r 1 ,r 2 ) 



i\ 2 P **(L> 



(ri,r 2 ). 



(10) 



The limits of the integral only cover the overlapping region. We notice here that if we use the Limber approximation Eq. JB2t . then this equation reduces 
to simpler form as higher harmonics at different radial distances n becomes uncorrelated. 



Cr(ki,k 2 ; ri ,r 2 ) = S 1D ( ri - r 2 ) 



2rf 



l 2 (l + l) 2 



ki k 2 Ii(ri,r 2 ). 



(11) 



The convergence power spectrum now can be computed using Eq.l[6]l. It is worth mentioning that this equation establishes a direct link of convergence 
power spectra and the underlying mass distributions. In later sections we will extend this result to higher-order multispectra. 

The signal-to-noise associated with various estimators from all-sky weak lensing surveys and the issues related to optimization will be dealt with 
in a separate work. In this paper we will focus mainly on development of statistics which can be employed to study gravity-induced non-Gaussianity 
using higher-order statistics. We will use the expressions for the Cis derived here for the construction of the optimal estimators for the bispectrum and 
trispectrum in the following sections. 



3 MODELLING GRAVITY-INDUCED NON-GAUSSIANITY 

Two point statistics are useful to constrain cosmological parameters. However the convergence power spectrum depends principally on a specific com- 
bination of cosmological parameters. Typical additional inputs in the form of external data sets such as the CMB, or higher-order statistics can lift the 
degeneracy. While the use of the 3D power spectrum already tightens the constrains further gain is anticipated with the help of non-Gaussianity studies 
in 3D. 

It is already known from earlier studies dTakada & Jainll2004T) that lensing tomography with the power spectrum and bispectrum can act as a probe 
of dark energy and mass power spectrum. The lensing bispectrum has a different dependence on the lensing weight function and the growth rate of 
perturbations. This is the main reason why bispectrum tomography can provide complementary constraints to the power spectrum. In fact it was found in 
previous studies that constraints from the bispectrum can be as tight as that from the power spectrum. 

We will model the non-Gaussianity using the hierarchical ansatz which is known to be a reasonable approximati on at small scales in the highly 
nonlinear regime. This approach has been used previously to model the statistics of the convergenc e and shear fields jMunshi & Jainll200ll :IValageas 
l2000l : lMunshi & Valageasl(2005l : IValageas. Barber. & Munshill2004l : IValageas, Munshi & Barbedl2005l) . 
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3.1 The Hierarchical Ansatz in the Highly Non-linear Regime 

We need a reliable technique for modelling weak lensing statistics beyond the power spectrum. On larger scales, where the density field is only weakly 
nonlinear, perturbative treatments are know n to be valid. For a s tatistical description of dark matter clustering in collapsed objects on small scales, the 
standard approach is to use the halo model (Cooray & Seth 2002). However, an alternative approach on small scales is to employ various ansatze which 
trace their origin to field theoretic techniques used to probe gravitational clustering. For our work, we will use a hierarchical ansatz where the higher- 
order correlation functions are constructed from the two-point correlation functions. Assuming a tree model for the matter correlation hierarchy (typically 
used in the highly non-linear regi me) one can write the m ost general case, the N point correlation function, £jv(ri, . . . , r„) as a product of two-point 
correlation functions £(|r-j — r 3 -|) dBemardeau et alir2002h . Equivalently in the Fourier domain the multispectra can be written as products of the matter 
power spectrum Pj(fci). The temporal dependence is implicit here. 

(JV-l) 

dN(rx,...,r a ) = (6(r 1 )...6(r n )) c = £ Q N , a ^ [] (12) 

a, N— trees labcllings cdgcs(i,j) 

It is however very interesting to note that a similar hierarchy develops in the quasi-linear regime at tree-level in the limiting case of vanishing variance. 
However the hierarchical amplitudes become shape-dependent in such a case. Nevertheless there are indications from numerical sim ulations that these 
amplitudes become configur ation-independent again as has been shown by high resolution studies for the lowest order case Q 3 = Q JScoccimarro et all 
1998; Bernar deau et a 1 2002). In the Fourier space however such an ansatz means that the entire hierarchy of the multi-spectra can be written in terms of 
sums of products of power spectra with diferent amplitudes Qn.c, etc. , e.g. in the lowest order we can write: 



<<5(ki)<5(k 2 )> c = (27vf8 3D (k 1 + k 2 )P(fci) (13) 
(<5(k 1 )<5(k 2 )5(k 3 )) c = (27r) 3 S 3 D(ki + k 2 + ks)B(ki,ka,ks) (14) 
<5(ki) ■ ■ • <5(k 4 )) c = (27r) 3 5 3 c(ki + k 2 + k 3 + k4)T(ki, k 2 , k 3 , k 4 ). (15) 

The subscript c here represents the connected part of the spectra. The Dirac delta functions Sgo ensure the conservation of momentum at each vertex 
representing the multispectrum. 



B«(ki,ka, k 3 )£ k . =0 = Q 3 [P(k 1 )P(k 2 ) + P{kx)P(k 3 ) + P(k 2 )P(k 3 )] (16) 
r 5 (ki,k 2 ,k 3 , k 4 )^ ki=0 = R a [P(k 1 )P(k 2 )P(k 3 ) + cyc.perm.] + I^[P(ki)P(\ki + kaDPflkj + k 2 + k 3 |) + cyc.perm.]. (17) 

Different hierarchical models differ in the way numerical values are allotted to various amplitudes. Bernardeau & Schaeffer dl992l) considered "snake", 
"hybrid" and "star" diagrams with differing amplitudes at various order. A new "star" appears at each order, higher-order "snake s" or "hybrid" diagrams 
are built from lower-order "star" diagrams. In models where we only have only star diagrams dValageas. Barber. & Munshll2004l) the expressions for the 
trispectrum takes the following form: 

T(ki, k 2 , k 3 , k 4 )^ k . =Q = Q i [P{k 1 )P{k 2 )P{k s ) + cyc.perm.]. (18) 

Following Valag eas. Barber. & M unshi < |2004» we will call these models "stellar models". Indeed it is also possible to use perturbative calculations which 
are however valid only at large scales. While we still do not have an exact description of the non-linear clustering of a self-gravitating medium in a 
cosmological scenario, theses approaches do capture some of the salient features of gravitati onal clustering in the highly non-li near regime and have been 
tested extensively aganinst numerical simulation in 2D statist ics of convergence of shea r ( I Valageas. Barber. & M unshi 2004). These models were also 
used in modelling of the covariance of lower-order cumulants (Munshi & Valageas 20051). 



4 THEORETICAL MODEL LING OF 3D CONVERGENCE BISPECTRUM FOR 3D WEAK LENSING SURVEYS 

Previous studies of the bispectrum involving weak lensing include work on projection (2D) as well as tomography. In most of these studies the prediction 
of the bispectrum is tied to a specific assumption regarding the growth of instability in the underlying density distribution. The main motivation for most 
of these studies was to put tighter constraints on the dark energy equation of state using weak lensing surveys by lifting the degeneracy involved in power 
spectrum analysis alone. The three-point correlation function in real space (or equivalently its harmonic transform the bispectrum) encodes information 
about the departure from Gaussianity, and this departure can be induced by non-linear gravity or by non-Gaussian initial conditions. We will focus on the 
gravity-induced bispectrum here and plan to present a complete treatment of the effect of initial non-Gaussianity on 3D weak lensing statistics elsewhere. 



4.1 Linking the Density Bispectrum in Various Representations: <5(k), 5i m (r) and 8i m {k) 

The convergence bispectrum in 3D will depend on modelling of underlying density bispec trum, We start by quoting the relation of the 3D density 
bispectrum expressed in Cartesian coordinate and in the harmonic space, we refer the reader to lCastro et all 12005) for detailed derivation of the following 
equation: 

5 lm {k;r) = -^ki l I dn k S{k;r)Y lm (n k ). (19) 

V27T ./ 
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Table 1. Notations in different bases for the bi- and trispectrum for the underlying mass distribution 8. The respective multispectra for the convergence field k will be denoted 
by corresponding calligraphic symbols Bi-ii 2 i 3 and 7jj 4 . 



Function 


3pt Correlation 


4pt Correlation 


Basis/Space 


S(r) 


£3 (ri,r 2 ,r 3 ) 


£ 4 (n, • • ■ ,r 4 ) 


Real Space 




B mixod (riir2ir3) 


T>^1^2 / r . . \mixcd 




Si m {k;r) 








<5(k;r) 


S rect (ki,k 2 ,k 3 ) 


T rect (ki,k2,k3,k 4 ) 


e ikx 



Using this definition, we can link the bispectrum defined in Fourier space with the one in the harmonic domain: 

(Si imi {ki;ri)8i 2m . 2 (k 2 ; r 2 )Si 3rri3 {k 3 ; r 3 )) c = 

(y=) {kik 2 k 3 }i ll+l2+l3 J dCl kl Y hrni (Cl kl ) J dh k2 Y hm2 (Cl h2 ) J dU k3 Y l3m3 (fl k3 )(5^r,ri)6(k 2 ;r 2 )S^3;r s ))^ (20) 
Let us introduce the following notation for the bispectrum associated with the density field: 

(<5(ki,ri)<5(k 2 ,r 2 )5(k3,?-3)} = B(ki, k 2 , k 3 ; n, r 2 , r 3 )S 3D (k 1 + k 2 + k 3 ). (21) 
Expanding the Dirac delta function 5 3 D(k± + k 2 + k 3 ) and using Rayleigh's expansion of the exponentials: 

S 3D (kl+k 2 + k 3 ) = ^3 J e >("l+k2+k 3 ).r d 3 r 



E 

Li, Mi 



d i vi Ll+L2+L3 jl 1 (kir )]l 2 (k 2 r)jL 3 (k 3 r)Yz /1 M 1 (f2fc x )Yl 2 m 2 (^fc 2 )Yl 3 m 3 (^fc 3 )Yl 1 m 1 (Q)Yl 2 m 2 (Q)Yl 3 m 3 (n). (22) 



Next, we use the orthogonality property of the spherical harmonics, Eq. dB5b to carry out the integrals to simplify the expression. We have introduced 
the notation <i 3 k = k 2 dkdCl^ = k 2 sin8kdkd8kd(f>k. This allows us to write the bispectrum in a spherical harmonics representation to its Fourier 
counterpart. Spherical coordinates are natural choice for various reasons. The line of sight integration can be treated quite separately with sky coverage 
issues. As we will see the partial sky coverage issues can also be dealt with a natural way in harmonic expansions. It is also important to notice that 
(radial) errors due to photometric redshift can also be incorporated naturally (see Appendix [A] for more details). 

{8i imi (ki;ri)) ...5i 3m3 (k 3 ;r 3 )) = y-^=) GTJ^ 3 J r2dr Jh( k i r )ji2( k 2r)ji 3 (k 3 r)B{ki] ri ). (23) 

The directional dependence through the azimuthal quantum number m is encapsulated through the Gaunt integral G and is defined by the following 
expressions (we also introduce the quantity Ii^^ which we will find useful later): 

GTjzr 3 = / dCiY hmi (Ci)Y h rn 2 (m 3m3 m = [ l i 1 2 „ i?„ Kw*; <^> 



mi m 2 m 3 



^ (2Z 4 + 1)(2Z 2 + 1)(2Z 3 + 1) I h h h 
lllhl3 ~ V 4tt I 



(25) 



Here the matrices are the 3J symbols, which are nonzero only if the triplets of harmonics (h, l 2 , l 3 ) satisfy the triangle equality, including the condition 
that the sum h + l 2 + l 3 is even which ensures the parity invariance of the bispectrum. We will also need the shorthand notation li 1 z 2 i 3 in our following 
derivations. The rotationally invariant bispectrum Bi 1 i 2 i 3 can now be written in terms of -B™^ 2 '™ 3 as: 



Br^r s ^r h = (^ t 2 t) B ^ {k ^ ^^^= E X ™ 3 ) B ^ m3(fc;r,) 



sph (26) 



m l ! m 2 j m 3 



We will also need the reduced bispectrum commonly used in the literature which has direct correspondence to the flat-sky bispectrum. In terms of the 
bispectrum Bi 1 i 2 i 3 we can define bi t i 2 i 3 — ^ 1 ! 2 i 3 -Bi 1 i 2 i 3 . Finally we can write the general correspondence between the spherical harmonics representa- 
tion of the angular bispectrum Bi 1 i 2 i 3 (fa;r-j) = Bi 1 i 2 ; 3 (fci, k 2 , k 3 ;n,r 2 ,r 3 ), and its Fourier representation Bi 1 i 2 i 3 (ki; n) = B(ki, k 2 ,k 3 ;ri,r 2 ,r 3 ). 
We will suppress the explicit display of the radial coordinates to simplify notations. We need to keep in mind in the Fourier representations the radial 
coordinates simply denote the cosmic epoch. 
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{fcife 2 fc 3 } Ji 1 i 2 i 3 (k l ;r i )I h i 2 i 3 



sph 



(fc<; ri) 



{kik 2 k 3 } Ji 1 i 2 i 3 (k t ;r i ). 



The expression Ji 1 i 2 i 3 (ki; fi) encapsulates the dependence on ki and n with (i = 1, 2, 3) 

Jiii-2i 3 {ki\ri) = / r 2 drB rect (k 1 ,k 2 ,k 3 ;r 1 ,r 2 ,r 3 )j h (k 1 r)ji 2 (k 2 r)ji 3 (k 3 r). 



(27) 



(28) 



In certain applications it is also interesting to work in a basis where the harmonic decomposition is carried out only on the surface of the sky, retaining 
the radial dependence in configuration space. In such circumstances, the following transformations are useful in relating the bispectrum expressed in this 
mixed coordinate with bispectrum in full spherical coordinate. 

&lm(r) = \J^ J kdkji(kr)8i m (k); (ri, r 2 , r 3 ) = (J^j J dk 1 k 1 j h (k 1 r 1 ) . . . J dk 3 k 3 ji 3 (k 3 r 3 )B B l ^ 2h (k 1 , k 2 , k 3 ;ri,r 2 ,r 3 ). (29) 

The inverse transformation from the basis <5; m (r) to <5; m (fc) and the related change in bispectrum are as follows: 

SUk) = sfl J r 2 drkji(kr)S hn (r); B*^ h (fei, fe, te; n, r 2 , r 3 ) = (J^j J drirfkij^ (feir) . . . J dr 3 rlk 3 ji 3 (k 3 r 3 )Bi^^ (n, r 2 ,r 3 )(30) 
If we put the expression Eq.J27t into the equation Eq.<|29t then we can finally write: 



B hhh( ri ' r2 ' r 3) = (J^) I hhh I dk±k 2 j h (k 1 r 1 ) . . . I dk 3 k 3 j h (k 3 r 3 ) I r 2 drj h (k 1 r)...ji 3 (k 3 r)B rcct (k 1 ,k 2 ,k 3 ;r 1 ,r 2 ,r 3 ) 



(3D 



These relations, e specially Eq.J27t will be useful in linking the 3D bispectrum to the convergence bispectrum. This is an extension of earlier results in 
ICastro et al (2005) for the power spectrum, where they showed Ci (k; r\,r 2 ) = P(k; n, r 2 ). The result at bispectrum level is more involved. In the next 
section we will use some approximations to simplify these expressions. 

4.2 Linking the Convergence Bispectrum to the Underlying Matter Bispectrum 

To make contact with the observables we use the fact that the projected convergence (which is related to the lensing potential) can be related directly to 
the 3D density field. We will start by linking the 3D convergence bispectrum B and the 3D density bispectrum expressed in harmonic coordinates. In the 
next section we will express the bispectrum in spherical coordinate in terms of the bispectrum in rectangular coordinates and use some well-motivated 
approximations to simplify the results. Using Eq.([5} we can write the following expression: 



dr 2 r 2 j h (fe 2 r 2 ) 



dk[ 

"fcT 



dr' 2 



r 2 - r 2 



dnnji 2 (k[r[) 

.1,1 

dk' 3 



ri dr\ 



T\ - r 1 



dr 3 r 3 ji 2 (k 3 r 3 ) 

Jo 



dr' 3 



r 3 



(32) 



The bispectrum Bi 1 i 2 i 3 (ki; ri) sph is now expressed in terms of the bispectrum B^^ ls (ki; ri). This relation mixes modes only in radial directions. On 
the surface of the sky there is no mixing of angular harmonics corresponding to various I values. While expressing the density harmonics in terms of the 
3D potential harmonics, we pick up additional scale factor o(r^) and ki dependence in the denominator (see Eq.([4j for more on notational details). 

We have so far ignored the presence of noise. Indeed because of the limited number of galaxies available it may not be possible to probe individual 
modes of the bispectrum at high signal-to-noise ratio. In later sections we will be able to address issues related to optimum combinations of individual 
modes which may be better suited for observational studies. 



4.3 Specific Models for underlying bispectrum and Limber's Approximation to the Exact Analysis 

The above analysis relates the underlying 3D bispectrum to the 3D convergence bispectrum for the most general case. However numerical computation 
involving such multiple integrals can be prohibitive. To make further progress we will use specific models of gravity-induced bispectrum to simplify the 
calculations. We will also use the Limber approximation to simplify our results. The Limber approximation is known to be a very good approximation for 
smaller angular scales or high I. We would like to stress however that although the results presents here are for a specific models for hierarchical clustering 
it is nevertheless possible to extend the results of our analysis to other models too. Assuming the specific form of hierarchical ansatz,introduced before, 
we can have: 



B roct (ki,k2,k3;ri,r2,r 3 ) = Q 3 [P(ki,n)P(k 2 ,r 2 ) + cyc.perm.]. 

Using this notation for the function Ji 1 i 2 i 3 (n, r 2 , r 3 ) we introduced in Eq.J28t takes the following form: 

Jhi 2 i 3 (ri,r 2 ,r 3 ) = Q 3 I h i 2h / r 2 dr j h (k 1 r)j h (k 2 r)ji 3 (k 3 r) [P(kv, ri)P(fc 2 ; r 2 ) + cyc.perm.] . 



(33) 



(34) 
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We use the extended Limber approximation (see Eo. dB2t ) to simplify the integrals involving k' i . The delta functions simplify the resulting r' integrations, 
Finally the observable convergence bispectrum can be written in terms of directly the density bispectrum as follows: 



B, 



Zi 1 i 2 i a (r 1 ,r 2 ,r :i ) 



dr 1 r 1 j h (k 1 r 1 ) 



dr 



dr 2 r 2 j h (k 2 r 2 
(n - r 



dr 3 r 3 ji 3 (k 3 r 3 )l h i 2h (n, r 2 , r 3 ) 



a(r)r 3 



<T2 



a(r)r 3 



a(r)r 3 



Ml 



r J P I — ; r J + eye. perm. 



(35) 



The integral here extends to the overlapping region i.e. r m in ~ min(ri, r 2 ,r 3 ), and the final result is not specific to the assumed non-Gaussianity, but 
assumes the hierarchical ansatz. For an arbitrary bispectrum the result can be expressed by a suitable change in Ii 1 ; 2 ; 3 (ri , r 2 , r 3 ) : 



Zhi 2 i 3 (ri,r 2 ,r 3 ) = 



l hhh 



dr 



(n - r) 



a(r)r 3 



(r 2 - r) 



a(r)r 3 



(r 3 - r) 



a(r)r 3 



B[—,—, —;r,r,\ 
V r r r 



(36) 



For computation of the bispectrum in scenarios with a specific model for the primordial non-Gaussianity we will have to replace the kernel that 
appears in the expression for Ii 1 i 2 ; 3 (ri, r 2 ,r 3 ) and similar results will follow. In particular we can replace the gravity-induced bispectrum with models 
of the primordial bispectrum, e.g. B loc or B equl , to compute the related bispectrum for convergence k. 



5 THEORETICAL MODELLING OF THE CONVERGENCE TRISPECTRUM FOR 3D WEAK LENSING SURVEYS 

As before we start by linking the trispectrum in spherical coordinates with the spatial trispectrum in rectangular coordinates. The procedure we will 
follow will be very similar to what we have done for the case of the bispectrum. We start by introducing the trispectrum in the Cartesian coordinate 
{5 (hi; 7"i) . . . S(k4\ r±)) — T Icct (k{\ ri)5 3 rj(k-i + k 2 + k 3 + k4) and in radial and polar coordinates as: 

(^ 1 m 1 (fci;r-i)5i 2m3 (fe;r 2 )<5 !3m3 (fc3;r3)5 i4m4 (fc4;r4)} c = J^(-1) M f ^ ^ M jf ^ ^ j T^fc (L, h; n) Bph . (37) 

LM ^ ' \ ' 

The vectors h, h, h, U represents the sides of a quadrilateral and L is the length of the diagonal. The matrices as before are the Wigner 3J symbols. The 
symbols are only non-zero when they satisfy several conditions; which are |/i — Ja| < L < h+l 2 , \ls — Ia\ < L < l 3 + h; h + 1 2 + L= even, l 3 + U + L 
= even and mi + m 2 = M as well as m 3 + rrii = —M. 

In our notation for the trispectrum, T^* [ki, ri\L), the indices (ki,ri) = (ki, ri, . . . , ki, r^) encodes their dependence on various Fourier modes 
of the density harmonics in the radial direction, used in their construction. No summation will assumed over these variable unless explicitly specified. We 
need also to subtract the Gaussian or the disconnected part from the estimated trispectrum to compute the connected part of the trispectrum, denoted by 
the subscript ( ) c in ensemble averaging. By expanding the Dirac delta function in spherical harmonics and going through the same algebra as above we 
can finally express [L, h; n] 3ph in terms of T A (h; n) Icct : 

\^(^-\\ M n'l'2 L n l 3 l i L I ^A~.A, A, (U.~..\T l l l 2(T,l.. ^ s P h 



(5 limi (ki, ri ) . ..S Um4 (k 4 ,r 4 )) c = }^(-l) M G^ 2M G^ 3 %_ M / rUrj h (k iri ) . . . ^(k^T^L; h, n) Bpi " (38) 

LM 



Next we express the four-point correlation function in terms of the trispectrum T* 1 ^ [L, ki\ ri\. Finally using the orthogonality properties of the 3J 
functions, we can finally connect the two representations. It involves the functions /; 1 ; 2 i 3 we have introduced before. The prefactor involving ki is an 
artifact of the normalization which we have adopted. 



4 



T^[L;ki,n] spb = ( -j= ) {k!k 2 k 3 ki} ^I lll2L I hliL J 4 (h;n); J 4 (fe;n)s I r 2 drT4{h;r t y cct j l (k 1 r)...j l (k i r) (39) 



" i ; ~ ' ; ■ f. ,. '-I* _ / _ I / z-. i — u~u ,\ r, ■ /. - ,. i- i.n i _ / 

L 

In our derivation we have used the following identity to simplify the integration involving four spherical harmonics: 

dQY limi (n)Y l2m2 (Q)Y hm3 (n)Y l4mi (^) = ^(-l) M G^ 2M Gg^- M (40) 

LM 

We will also add the following expressions for the sake of completeness. As before we will relate the trispectrum defined from the harmonics 5i m (k; r) 
i.e. Tjgj 2 (L, ki\ri) with T t ^ t 3 (L; ri) which is defined from the harmonics Si m (k). 

T\}Jl (L; ki,nr h = (I ) 2 J rUnk^ fan)... J rjdr^ju (k^T^ (L, r^™^ (41) 
The inverse relation which relates T ; ^ (L, ri) with (L; ki,ri) is given by following expression: 

Tllll (L, r t ) mixcd = ( I) 2 J fcirffciji (feir) ...J ktdluji (^r)T^ (L; k^n)^ (42) 

The relation of the full 3D trispectrum in spherical coordinates and its Fourier decomposition, which generalizes our previous results for the bispectrum, 
is written as follows: 
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T ilil (L,ri) m,7ced = (J^j y^Jhi 2 Lli 3 UL J dkikij h (kiri) . . . J fcldfc 4 j i4 (fc 4 r 4 ) J r 2 drji 1 {k 1 r 1 )...ji 4 (k 4 r i )Tl^(L,k i ,ny cct . (43) 

Along with the bispectrum expression this generalizes the previously-obtained relationship at the level of the power spectrum in ICastro~et al (2005). 
Clearly for practical purposes we will need to devize an approximation to the exact result. We will use the Limber approximation to approximate the 
spherical Bessel functions. Note that numerical evaluation of trispectra is considerably more involved than the bispectrum. It is also important to note that 
as we climb upwards in the hierarchy realistically it gets more difficult to extract signals from observational data because of the presence of noise. 



5.1 Linking the convergence trispectrum with the underlying matter trispectrum 

Finally the observable trispectrum T for the convergence ft (defined through an equivalent expression as in Eq.ll37ll) can be expressed in terms of the 
underlying trispectrum of the mass distribution: 

Ti^i* {h; ri) sph = kik 2 k 3 ki J jj- J dr 2 r 2 ji 2 (k' 1 r' 1 ) J 



TlllUK-r'^. (44) 

The mode-mixing in spherical coordinates happens only in the radial direction. It is expected that the estimation of the trispectrum from a realistic sky will 
be noise-dominated in the near future. This means estimation will be difficult for individual modes. We will develop methods to compress the information 
content in individual modes in an optimal way elsewhere. The trispectrum is dominated by the noise from galaxy intrinsic ellipticity as well as shot noise 
from the Poissonian nature of the galaxy distribution. To determine this we need to take into fact that a contribution to the trispectrum not only comes 
from the non-Gaussian signal but also from disconnected Gaussian terms too. 




r dk' A r°° . , , r 

Jo K t Jo Jo 



5.2 Specific Forms for Underlying Matter Trispectrum and the Limber Approximation 

We will derive the result quoted above for the case of the hierarchical ansatz with a "stellar" approximation we make further use of the extended Limber 
approximation to simplify. The hierarchical ansatz as well as the Limber approximation are both valid at the small angular scale, which justifies their joint 
use to simplify the results. The result takes the following form: 



dnriji 1 (kin) 



drAVij^^kir^Xi^isii (n, r 2 , r 3 , r 4 ), 



li 1 i 2 i a i i (r 1 ,r 2 ,r 3 ,r i ) 



E 

L 



hii 2 Lh 3 i i L 



1 dr 


(n - r) 




(r 4 - r) 


Jo 


a(r)r 3 




a(r)r 3 



(45) 



These results are extensions of analogous relations obtained for the bispectrum. We will introduce contributions from star topology under the stellar 
approximation (for other hierarchical ansatze see e.g. lSzapudi & Szalavl 11993LI1997I) which assumes that the amplitudes associated with all topologies 
are the same). 

We are only concerned with the connected part of the trispectrum here. Next we use the hierarchical ansatz to model the four-point correlation 
function. We will only use the contribution from the diagram with "star" topology in this section. The trispectrum expressed in Fourier domain, is simply 
the Fourier transformation of the four-point correlation function. The trispectrum as outlined before in hierarchical approximation can be written as a 
product of three power spectra: 



(<J(ki, n)5(k 2 , r 2 )<5(k 3 , r 3 )c5(k 4 , r 4 )) c 



— Rq 

+R b 



d 3 kP{k 1 )P{k 3 )P(k)S 3D (k 1 + k 2 - k)<5 3D (k 3 + k 4 + k) + cyc.perm. 
d 3 kP(fci)P(fe).P(fe 3 )<5 3D ((ki +k 2 -k)5 3D (k 3 + k 4 + k) + cyc.perm 



(46) 



In general the hierarchical amplitudes R a (associated with the snake topology) and Rb (associated with star topology) will have different amplitudes. 
There are 12 terms with snake topology and 4 terms with star topology which are represented by the "cyc.perm". Various hierarchical mod els differ in 
the way they ascribe values to various amplitudes. It is possible also to employ Hyper Extended Perturbation Theory dScoccimarro et alfl 998) to compute 
these amplitudes. For our purpose we will assume: 



(<5(ki; n)5(ka; r 2 )5(k 3 ; r 3 )<5(k 4 ; r 4 )) c = Q 4 



d 3 kP(fci)P(fc 3 )P(fe)<5 3D (ki + k 2 - k)<5 3D (k 3 + k 4 + k) + cyc.perm. 



(47) 



In this case, the analysis is essentially the same as that of the bispectrum and it simplifies the results considerably. The stellar approximation consists 
of approximating the four-point correlation only with stellar diagrams. This model has been checked in considerable detail in 2D in previous work 
tearber. Munshi & Valageasl 2004: Mu nshi. Valageas & Barbed 1 20041 : IValageas. Munshi & Barber! 120051) . We will assume a "stellar" model from this 
point onward. However the method outlined can also be generalised to take into account the "snake" diagrams. The cyclic permutations for the stellar 
model now represents all 16 diagrams. 
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We start by expanding the Dirac delta functions 630 (k) using two dummy positional variables x and y. Next following the same procedure as we 
have followed for the case of the bispectrum we can express the star part of the trispectrum in spherical coordinates. This will next be needed for the 
derivation of the convergence trispectrum. 



{8i imi (ki;ri)...5i imA {kA\r4))o= / drmj h (fan) . . . / dr 4 r 4 j i4 (fc 4 r 4 ) J iii 2 i 3 i 4 ( ri ' r ' 2 > r3 ' r4 )' 

Jo Jo 

JlllsU^'^rs,^) = Q 4 (^) 4 fc 1 fc 2 fc3fc4^G^ M G™- M f r 2 drj h (k 1 r)...j l4 (k 4 r){P(k 1 )P(k 2 )P(k,) + cyc.perm.}, (48) 



where the definition for the kernel J^l 2 i 3 i i (n, r 2, f 3, r 4 ) is similar to its counterpart we introduced for the bispectrum. In our notation, A l ^\* denotes 
the star contribution with the corresponding amplitude R a . It is now possible to express the star contribution to the 3D convergence trispectrum using the 
following relation: 



(fa;n) . . . Ki 4 m 4 (fa; r 4 ))c 



dr 1 nji 1 (fan) . 



drir4ji 4 ,(k4r i )li 1 i 2 i 3 i i (r 1 ,r 2 , r 3 , r 4 ) 



T ll i 2h i i {r 1 ,r 2 ,r 3 ,r A ) = C 



( 4 A) 

\TVC 2 J 



/4fa\ r7rl 
Vttc 2 / L 2 J 



E 

L 



1 dr 


(n - r) 




(r 4 — r) 


'0 


o(r)r 3 




a(r)r 3 



x(5 iimi (fa; n) ■ • • <5i 4 m 4 (fc4; r 4 )) st0 r- 



(49) 



We simplify the expression using the Limber approximation as before. It effectively replaces the wavenumber fc;s with the corresponding U/r while 
reducing the dimensionality of the integrals. 



(ti imi (fa; ri) . . . K i4m4 (fc 4 ;r 4 )) c = / dnnjix (fan) / dr 2 r 2 ji 2 (k 2 r 2 ) / drsrsjia^sni^iaia^i.^rs) 

Jo Jo Jo 



Ii 1 i 2 i 3 {ri,r 2y r :i ) = 



L 



dr 



(n - r) 
a(r)r 3 



(r 4 - r) 



a(r)r 3 



{ p (t"' 7- ) p (^' r ) p (^' r ) + C y°-P erm -} • ( 5 °) 



Here the upper limit of integration along the radial direction is rw n = min(ri, r 2 , rs, r 4 ). In our analysis above we have taken the hierarchical ansatz 
as an example, but it is quite general and the expression for the density trispectrum only affects the expression for Jj*j j ; 4 (ri, r 2 , rs, r 4 ). Perturbative 
results are in general more complicated to deal with because of configuration angle-dependence, but will result in similar signal-to-noise ratio. In stellar 
models the star topologies at various orders carry all the weights, diagrams with snake topologies are ignored and arbitrary order in correlation functions 
are simply expressed in terms of the star contributions at that order. This approximation, as we will see, can simplify the calculations immensely. 

We have concentrated here in modelling of the trispectrum and stressed its importance in confirmation of detection of non-Gaussianity determined 
using the bispectrum. The analytical modelling of the trispectrum is also important in itself for calculation of the error-covariance of the power spectrum. 



6 CONVERGENCE SKEW-SPECTRUM 

Because of signal-to-noise issues it is difficult to study the bispectrum for all possible configuration of triplets. The skewness compresses all the infor- 
mation contents of the bispectrum into a single number. Such aggressive data compression may be elegant but it fails to d istinguish various contrib utions 
which might have different shape dependence. These issues have been discussed extensively in recent literature (see dMunshi & He avens 2009]) and 
references therein for detailed discussion of related issues). The cumulant correlators which were introduced in the literature are the two-point objects 
and are well studied in the case of galaxy surveys. These are muti-point statistics collapsed to two-point objects. In harmonic space they correspond to 
power spectra ass ociated with mu ltispectra of various orders. At the lowest order there is only one power spectrum (coined the skew-spectrum) related 
to the bispectrum. Coorav (2005) had earlier introd uced the unopt i mised versions of power spectra associated with bispectra and used them to study 
lensing-secondary cross correlation. Later studies bv lCooravl ( l2006l) : ICoorav. Li & Melchiorril d2003h used power spectra associated with the bispectrum 
and trispectrum for redshifted 21cm studies. These power spectra retain some of the information of the multispectra that they are associated with, and the 
number increases with the order. We will generalize and use the idea of cumulant correlators here to study the bispectrum and trispectrum associated with 
the 3D convergence fieldQ. 

The squared convergence field K 2 (ri) is constructed at a radial distance n; its harmonic transform on the surface of the sky is denoted by «S(n) 
(we will be using the mixed representation throughout). Let us start by expressing the spherical harmonics transform nf^ (n) of the squared convergence 
field K 2 (£l, n) in terms of the spherical harmonics of the original convergence map «i m (n). 

<4%{ri) = J r^(n) K 2 (n,n)df2= ^2 KJimitriKmafri) J dAyi imi (A)y Iaraa (A)y,^(ft). (5i) 



4 Detailed modeling of a multispectra is not important for denning the associated power spectra 
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The above expression assumes all-sky coverage. In practice the surveys will cover only a fraction of the sky. The mask used in the survey w(ti) will affect 
the estimators that introduces a multiplicative bias which needs to be properly accounted for. If we denote the masked sky harmonics of the squared field 
by (n) then we can express them in terms of the original convergence harmonics of the all-sky and the harmonics of the mask: 



J w(fl)Yi m (ti)n 2 (£l,ri)d£l = X] X/ K hm 1 (ri)Ki 2m2 (ri)wi 3m3 J dflY hmi (Cl)Yi 2m2 (Cl)Yi 3m3 (Cl)Y* m (D.) (52) 

limi l'ira 2 l-^rn-i 

(I'm') 

Here the mixing matrix -R"; m ;' m ' [w] denotes the mixing of harmonics modes due to the presence of the sky mask whose harmonic transform is wi m - We 
will use the same mask at both radial distances, but the results can easily be generalized to two different masks. Using these harmonics we can now define 
the skew-spectrum (or the power spectrum related to the bispectrum), C\ ' , as in Muns hi & Heavens! (2009). The cut-sky version of the skew-spectrum 

~(2 1) {2) 

is denoted by C\ ' and is constructed from the cut-sky harmonics and ki m as follows: 

C\ %1) {vx,r^= ^l T ^Real{ K W*(r 1 )« !m (r 2 )}; (n, r 2 ) = £ Real (n)ft«m(ra) } . (54) 

m m 

The skew-spectrum C, (ri, r 2 ) probes directly the bispectrum Bi 1 i 2 i 3 . Though it does not encode the entire shape dependence for each triangular 
configuration it encodes more information compared to its one-point counterpart the skewness Sa{r\,r2), and has the ability to distinguish different 
contributions to non-Gaussianity, both primordial as well as gravity-induced. We will focus on gravity-induced non-Gaussianity here and issues related 
to primordial non-Gaussianity will be addressed elsewhere. 



„(2,i), s k / \ /(2Zi + l)(2i2 + 1) ( h h h\ 1 T „ , \ 

C ; (r 1 ,r 2 ) = ^B llll2 {r 1 ,r 1 ,r 2 )J I Q Q Q I = — y ^ I Ulh Bu lh (n, n, ra). (55) 



4tt(2/ + 1) V 1 21 + ' 

Here Bi 1 i 2 i 3 (n 1 r2,r3) is the angle-averaged bispectrum. It encodes information about the three-point correlation function in the harmonic domain, 

{r{)Ki 2 m, 2 {ri)Ki 3 m 3 (r 2 )) c : 

(K hmi {ri)K hm2 (r 2 )Ki 3m3 (r3)}c = B hh i 2 (n,r2,r 3 ) ^ ^ ^ ^ J . (56) 

While the bispectrum is invariant under the permutations of the angular harmonics (h,h,h) it is not invariant under the permutations of the radial 
distances r\, r<2, r 3 . In certain circumstances the reduced bispectrum bi x i 2 i 3 (ri, r 2 , r 3 ) is also useful to convey equivalent information: 



, (2Zi4-l)(2k4-l)(2i 3 + l) ( h h h , . , • 

B h i 2 i 3 {ri,r2,r 3 ) = \ n n n ^2*3 ( r i> r 2, r 3 ) = Ii 1 i 2 i 3 b h i 2 i 3 {n,r2,r 3 ). (57) 



4tt ^000 
For partial sky coverage one can obtain after tedious but straightforward algebra: 



r(2,i), , ST,n,> (21" + 1) ( I V ("V, px^„ , \ / (2ii + l)(2f 2 + l) / ii i 2 J' 



J_ r.a. i„ |2 



- 2/ + , ^/zVK"! < gp+I Z^ /l ''i , » B, W» , i,ri,f*) f = ^H'<r (n.ra)- (58) 

I hh ) 

wi is the power spectrum of the mask, which is completely general. It is easy to see that in the absence of any correlation between signal and noise, the 

(2 11 

estimator C\ ' is unbiased and no noise subtraction is needed as long as it is Gaussian. Using the definition of the coupling matrix M, introduced above, 
we express the pseudo-Cl (PCL) estimator as: 

Cf' 1) (r 1 ,r 2 ) = A/-, 1 4 2 ' 1) (r 1 ,r 2 ). (59) 

The covariance properties of such an estimator can be computed using similar techniques. If we collapse the two-point correlator to a one-point object, 
we can write the cross-skewness as 

&(n,ra) = £}(2Z + l)C? ll) (ri,r 2 ). (60) 
l 

In our analysis we have so far only used the statistics (k 2 ' (ri)K,(r 2 )) , which probes the bispectrum Bi 1 i 2 i 3 (ri,r\,r 2 ). It is however also possible to 
consider the analogous statistics (k(Ci, n)ft(f2, r 2 )n(Cl' , r 2 )) which probes Bi 1 i 2 i 3 (ri,r 2 ,r 2 ). There are similar other permutations which provide 
complementary information, with almost identical analysis. 
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7 KURT-SPECTRUM, OR THE POWER SPECTRUM ASSOCIATED WITH THE TRISPECTRUM 

The four-point correlation function, or its harmonic c ounterpart the trispectrum, has be en studied in the literature extensively. This contains the information 
about the non-Gaussianity beyond the lowest level dHull 19991 : lOkamoto & HvJl2002h . For the case of weak lensing studies clearly the gravity-induced 
non-Gaussianity is the main m otivation. Studies in trispectrum analysis have also been pur sued using various other probes e.g. using 21cm surveys 
iCoorav. ET & Melchiorri 2003) or more extensively in several CMB studies; see Bartolo et al ( 20041) for a review. However these studies typically probe 
the trispectrum induced by primordial non-Gaussianity. 

It is important to note however that at the level of four-point studies, the Gaussian fluctuations from the signal as well as the from the (Gaussian) 
noise too carry a non-zero (unconnected) trispectrum. This degrades the signal-to-noise for various estimators and clearly needs to be subtracted out 
before an unbiased comparison with the theoretical predictions can be made. 

It is obvious that detection of the trispectrum from noisy data is far more nontrivial than the estimation of the bispectrum. Previous studies have 
mainly concentrated on one-point estimators which collapse the data to a single number - known as the kurtosis. We extend studies involving kurtosis 
{k 4 (Q)} to its two-point counterparts: (ft 2 (f2, ri)K 2 (Cl' , 7-2)} and (n 3 (fl, ri)K,(£l', ^2))- In practice however we will consider the Fourier transforms of 

(3 1) (2 2) 

these objects which are the power spectra associated with the trispectra, C\ ' and C\ ' . Indeed the radial coordinates associated with two different 

(3 1) (2 2) 

fields being cross-correlated can be different and will be denoted 

We start by defining the all-sky harmonic transform Ky^(ri) for the convergence field K 3 (fi,ri) and cross-correlate it against K,i m (r 2 ). In the 
presence of a mask w(Cl) which we assume to be the same at both radial distances, the harmonic transforms of the cubic field (n) will depend also 
on the spherical transforms of the mask wi m too: 



K im( r l) = ^2 ^2 «!imi( r l) K !2m 2 (n)«! 3 m 3 (n) / dQ.Y hmi (Cl)Y hm2 {Cl)Y hm3 (A)F(^(fi) (61) 

ll mi ^2 m 2 

^ ( m( r i)= ^2 ^2 ^i m i( ri )^2m2( r i)^3m 3 (n)^ 4 m4 J dClYi imi (Q)Yi 2rri2 (Q)Y i3m3 (Cl)Yi 4m , 4 (0)l7 m (Q) 

limx l2m 2 l3 m 3 

= ^^ m ,' m '4 3 i'(ri)- (62) 

I'm' 

~/ o -|\ ""(2 2) (2 2) 

We will use these results to derive expressions for C\ ' (n, r 2 ) . The other cut-sky power spectra C\ ' (r 1 , r 2 ) and all-sky counterparts C\ ' (ri , r 2 ) 
are given by the following expressions: 

f (2,2), x _ 1 (2)» (2) £(2,2) _ 1 V- ~(2)» -(2) ™ 

m m 

These power spectra directly probe T^ 2 (/, n, r2) m,xcd . It compresses all the available information in quadruplet of modes specified by (Zi, Z2, £3, 14) to 

(22) (31) 

a power spectrum. The power spectra C\ ' (ri,r 2 ) and C, ' (ri,r 2 ) differ in the way they associate weights to various modes. 



.(2,2), rphhr, ^mixcd / (2? 1 + l)(2/ 2 +Ty / (2Z 3 + 1) (2Z 4 ( h h I \ ( h h I \ 

C, {r u r 2 )- 2^ H 3 U^ r ^) V 4»r(2l + l) V 47r(2? + l) I oilo J 

^1 ^2 7^3 

= (^TTF ^ IiLiiJusuTj^l,r u r 2 r* ed - (64) 

Here the reduced trispectrum T^ t * (rn L} mlxed j s defined in terms of («:; imi (ri)ni 2m2 (r 2 )ni 3Tn3 {rz)i^i i m l {ta))c as follows. We have added the radial 
distances r, associated with each spherical harmonic in the argument with L, which specifies the diagonal formed by the quadruplet of four quantum 
numbers U. 

(K hmi (ri)K l2m2 (r 2 )Ki 3m3 (r3)Ki imi (r4)) c = ^ T/^ 4 (r t ; L) mlxcd ( ^ ^ ^ J ( ^ ^ _^ . (65) 

For partial sky coverage we can express the cut-sky version of the estimator C^ 2,2 ' (n, r 2 ) in the following way. The resulting pseudo-CiS can then be 
expressed in terms of C,, (ri,r 2 ) through the mixing matrix Mu* : 



v rphh n xmw /(2Ii + l)(2i 2 + l) /(2i s + l)(2J4 + l) f Zi « 2 I' \ f h h I' 

X i I,l 4 C' "ll r 2j 



' zy V 4tt(2Z' + 1) V 4tt(2Z' + 1) ^ ^0 

1 r2 ,2 f 1 1 



2/ + 1 ^ 

VI" 



J "''"i^"! 2 { 2TTl E ^Vl^^^^ ;4 Lr/ 3 ^(;,r 1 ,r2) mixcd } = ^ M ii ,C i (2 ' 2 '(r 1 ,r 2 ). (66) 
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There are two cumulant correlators at four-point level as explained above. Following the discussion above we now focus on the other degenerate power 
spectra associated with the cumulant correlator k 3 (Q,)k(Ci'). This is of the same order as k 2 (Ci)k 2 (Ci') and contains information about trispectra as well. 
The compression of the information is done with different weighting for different modes: 

CP' 1 ^!,^) = ^ T ^Rea]{«^(r 1 )^(r 2 )} (n,r 2 ) = ^L. ^ Rea l {s^VO^fa)} . (67) 

m m 

We can now use the definition of the trispectra Tf 3 f (L; n, r 2 ) to express C,' 3 ' 1 ' (ri, r 2 ) in terms of the trispectra. The main difference with the previous 

spectrum Cp' 2 ' (r%, r 2 ) is that, it sums over all possible configuration of the quadrilateral keeping one of the sides fixed, whereas C[ 2,2 ^ (r%, r 2 ) keeps one 
of the diagonal fixed but sums over all possible configuration of the quadrilateral. 



Ms,i), v \- rjtibf. V nixcd / (2h + l)(2b + l) / (2L + l)(2i 3 + l) / h Z 2 L \ L h I \ fgse . 

C ' (ri ' ra) = ii 2^ £ T '. 1 « a ( i ^^) V 4^(2L + 1) V 4 ^ + !) 1° oj^O oj (68) 

= 27TT S ^TT /ili2L/L!3!T ^ 2(i ' ri ' r2)miXCd - (69) 

*2 £ 3 J L 

The partial sky coverage will mean that the measured power spectrum Ci 3,1 ' (n, r 2 ) is not the same as theoretical expectation, but is related as before by 

3 3,1) (ri,ra) = M u/ Clf- l) (r 1 ,r 2 ). (70) 

In fact it can shown that for arbitrary sky coverage with arbitrary mask the above analysis can be generalized to arbitrary order of correlation hierarchy. 
If we consider a correlation function at p + q order, for every possible combination of (p, q) we will have an associated power spectrum. Using the 
same expression for the mode mixing matrix, we can invert the observed Cf' 9 (ri, r 2 ) to Cf' q (ri, r 2 ). Hence for arbitrary mask with arbitrary weighting 
functions the deconvolved set of estimators can be written as: 

(71) 

The Gaussian components of the corresponding multispectra at that order need to be subtracted out, and can be written in terms of the C;. The noise 
contribution too is assumed Gaussian and hence only contributes to the unconnected components. As before we can collapse the two-point objects and 
reduce them to a one-point number, the cross-kurtosis, which will be a function of both radial distances ri,r 2 , 

K 4 { ri ,r2) = ^2(21 + l)Cj 3 ' 1) {r 1 ,r 2 ) = J^(2Z + l)C ; (2 ' 2) (n, r 2 ). (72) 
i i 

As we demonstrated with the cross-skewness, K4 (n, r 2 ) can be decomposed in Fourier modes in the radial direction and an associated power spectrum 
can be defined. 

The Gaussian contribution to the trispectrum can be written as: 



G\\\l (n, ra, rs, r 4 ; L) = ^ (2 Ii + 1)(2Z 3 + l)C h (n,r 2 )C h (r 3 , ri )5 L0 6 hh 8i 2h 

+(2L + l)(-l) h+l3+L 5 h i 3 5 hh C h (n, r 3 )C h (I a , r 4 ) + (2L + l)C h (ri,r 4 )C h (r 2 , r 3 )5i 1 iji 2h . (73) 

(3 1) (2 2) I I 

Next we can compute the Gaussian contributions to C\ ' and C\ ' following the same procedure as before just by replacing the trispectrum T t ^ with 

to realize that in computing the Gaussian contribution we will have to take into account both the signal and the noise C;s (assumed to be Gaussian), i.e, 

3f + Of 



its Gaussian counterpart (r»; L). Indeed we will have to keep in mind the ordering correct for various U and their r% counterparts. It is also important 
;omputin 

a = cf + a N 



G ! (3 ' 1) (n,r 2 ) = - ^2 2L 1 + 1 Iii1 2 lIli 3 iG\1\ 2 (L, n, r 2 ) mlxcd 

^1 &2^3 



Gi(2 ' 2)(ri ' r2)= (27TTF / »^ i »3U^an J r 2 r i * ed (74) 

For realistic surveys with a mask, results identical to Eq.d701l and Eg.ll66b will hold true for the Gaussian contributions. From the estimated Cj 3 (ri , r 2 ) 
and C*P' 2 ' (n , r 2 ) these contributions need to be subtracted out before comparing them against the theoretical expectations. 



8 OPTIMAL ESTIMATES 

The estimators introduced above are not optimal as they are not inverse-variance weighted. In this section we discuss the signal-to-noise of the estimators 

(21) (31) (2 2) 

introduced above, namely C\ ' at the level of bispectrum and C\ ' and C\ ' at the level of trispectrum. 
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For construction of t h e optimum estimates the harmonics a\ m needs to be inverse covariance weighted dSmith & Zal darriaga 2006) . We refer to 



iMunshi & Heavens! d2009l) : lMunshi et"al] d2009T) for a complete discussion which also the requires presence of a linear term {Creminelli et alj|2006h in 
the case of absence of spherical symmetry. The optimal estimates at the three-point level in the all-sky limit and in the presence of constant variance noise 
can be expressed as: 

5 ; (2 ' 1) (ri,r 2 ) = 2[ 1 - y^ j Bu 1 i 2 (r 1 ,r 1 ,r2)Bi h i 2 (r 1 ,ri,r 2 ) {Ci(r 1 ,r 1 )~ 1 C h {r 2 ,r 2 )~ 1 Ci 2 (r 3 ,r 3 )~ 1 + cyc.perm} . (75) 

We have denoted the estimators by". The other terms can be obtained by circular permutation. Different choices of ri,r 2 , r 3 can give us different skew 
spectra. In our discussion of the unoptimized version we have developed one specific example which corresponds to Si(r\, r\,r 2 ) and was denoted by 

(2 ll ('2 1) 

Cj ' \r\,r 2 ). The other choices that we can construct are Si(r\, r 2 ,r±) and Si(r\, r 2 ,r 2 ). It is important to note that C{ ' (r\,r 2 ) is not invariant 
under permutation of its indices. 

Similar results hold for the case of trispectrum. The optimal estimates for the case of the trispectrum, we can write for the full sky estimates with 
uniform noise: 



r(3,i) 



(n ' r2) = WTl ^ 2rU f ^ 2{L ' ,ri)T ^ 2(L ' ,n) { C ' 1 (^^i) _1 C ;2 (r 1 ,r 1 )- 1 C ;3 (r 1 ,r 1 )- 1 C ; (r 2 ,r 2 )- 1 +c y c.perm.} (76) 
f' 2) (n,r 2 ) = ^^ f ltll^n)Tllll{l-,n) {C h (n, nJ-'C^n, n^C^n, r 2 )" 1 C !4 (r 2 , r.y 1 +c y c.perm.} . (77) 



The one-point estimator with inverse covariance weighting is simply the sum over the free index of the respective two-point estimators and is expressed 
as follows: 

&(ri,ra) = ^(21 + l)i, (a,1) (r0; &{n,ra) =^(2l + l)ty 3,l \ri) = ^(2l + l)£f* ) (r i ). (78) 
l i i 

Depending on various choices to identify the quadruplet of the radial distances r» that are associated with each harmonic U we will have a different esti- 
mator which can provide complementary information. These estimators can help us to optimize survey depth and width for the study of non-Gaussianity 
at a given order. 



9 FLAT SKY TREATMENT 

As pointed out before for surveys with large opening angle, the all sky expressions developed so far involve the expansion in terms of the spherical 
harmonics and spherical Bessel functions. However for surveys which cover only a small patch of the sky most of the signal comes from higher harmonics. 
In such a situation, a natural choice would be to directly deal with a two-dimens ional Fourier expans ion suitable for flat space. This makes the analysis 
more straightforward. Our analysis here closely follows that of lHuJ d 19991) and lOkamoto & Hul d2002h . We are however required to take into account the 
additional radial coordinate in our analysis. We expand the 3D convergence field as before both in the line-of-sight direction as well as on the surface of 
the sky: 

k ) = \l~ J dr \\ J -^- K ( r h r ^) k ji( kr \\) exp(-il • rj_). (79) 

In this expansion 1 depicts a 2D angular wavenumber and k represents a wavenumber in the radial direction. It will also be advantageous in certain 
situations when a harmonic expansion is only performed on the surface of the sky, but the radial dependence is kept in configuration space. As before 
we have assumed in the above expansion that Universe is flat. Alternatively the eigenfunctions for the expansion needs to be suitably modified. The real 
space correlation functions and their Fourier counterparts are related by the following expression: 

/j2| ^ 2 j 
' ' ' ~2^~ r,ll ' ) ' ' ' K< " ln ' r ""^ c exp ^ ' r±1 H h ln ' rx ™)]- (80) 

The flat sky correlation hierarchy, which ensures translational symmetry of the 2D patch sky is expressed by the following equations. The r label is 
retained as we have not the performed the Fourier expansion in the radial direction. 

(K(h,r h )K(h,r ]h )) c = (2tt) 2 5 2D (1 1 + \ 2 )V{h, r B .) (81) 

( K (h, r y Xla, r ll2 ) K (h, r u )) c = {2-k) 2 8 2D (\ 1 + 1 2 + l 3 )B 3 (/i, h, h; r„ 4 ) (82) 

(K(h,r lh )K(h, r ]h Mh, r\\ 3 MU, r u )) c = (2tt) 2 5 2D (1 1 + 1 2 + 1 3 + l 4 )7i(/i, h, h, k; r h ). (83) 

The labels Ti which appears as the arguments for multispectra denote all the radial coordinates which are involved in their definition, e.g. rj = ri, . . . , r 3 
for the bispectrum. The treatment for trispectra is more complicated as it also gets disconnected Gaussian contributions (which are non-zero even in the 
absence of any non-Gaussianity) and the contribution from the reduced segment, discussed above, which carries all the information about non-Gaussianity 
at the level of fourpoint. 
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(K(Ii,r|| 1 )/s(l2 ) r|| 2 )K(l3,r|| 3 )/s(l4,r|| 4 ))G = (27t) 2 c5 2 d(1i + 1 2 )(27t) 2 5 2 ,d(13 + U)V(h, r\\ x )V{h, r\Q 

\2< 



+(27r)^«5 2D (li + l3)(27r)^ 2D (l 2 + U)T{h,r^)T{h,r {{ ) + (27r)^ 2D (li + h)(2nyS 2D {h + h)V(h,r lh )V{k,r u ) 
(K(li,n)K(l2,r 2 )K(l s ,r u )K(h,r u )) c = (2ir) 2 6 2 D(li +h +h + k)Mh,h,l 3 ,U,L;r h ). (84) 
Following the discussion in IhH d 1999b and lOkamoto & Hul (2002) we write the trispectra TJ ; in terms of the reduced trispectra as follows: 

T}£ = n\i\i(i 12 )+n\i\i{i 13 ) + n\i\-(h 4 ). (85) 

We have used the notation I12 = li + I2. We will quote the results obtained in lHul jl999h which relate the multi-spectra defined in the full-sky analysis 
to their flat-sky counterpart. 

C,(k;r) =V(k;r); B hl2h (n) = I llhh B(h,l 2 ,l 3 ; ri ); ^(L,n) = WW^^An). (86) 

While the results for the power spectrum and bispectrum are straightforward the reduced trispectrum is more involved as it depends on the choice of 
the diagonal of the quadrilateral constructed out of the vectors Us. The implemenation of the momentum conservation is imposed by the translational 
symmetry is built in the definition is 

[6 2D {h +h + l)S 2D {h + U - l)ftj££(J) + 4uj(1i + la + l)S 2D (h + U - (0 + 8 2D {h + 1 4 + l)<5 2 r.(l 2 + 1 3 - I)^jj(i)] d\ (87) 



/ 



The flat patch wave numbers (U) are used within the parentheses which appear on the r.h.s. of the equations, whereas their all-sky versions appear 
in the l.h.s. without the parentheses. 

The radial independence in our calculation can also be displayed by doing a Fourier transform in the radial direction. The transformations from 
all-sky to the Fourier representation are given by the following expressions. We use the same notations S3 or 71 in both representations. 

B 3 (li,l 2 ,l 3 ;n) = I J - ) dr^jiikir^) ... / dr\\ a ji(k 3 r\ ]d )Bz{h,h,h;ki). (88) 



A similar expression holds for the trispectrum: 

Ti{h,...,h\ri)= I J- ] dr ]h ji(kir tl ) ... / dr^ 3 ji(k 3 r\\ 3 )T4(li, ■ . ■ , U; h). (89) 



For a flat sky we can work with various representations for the multispectra as before. The real space variables r = (rii , rx) and their Fourier 
representations using variables k = (k, 1) are both useful. 

(js(fci,li)«i(fca, la)) = (2tt) 3 <5 2D (1i + h)5 1D (ki + k 2 )P(ki,h) (90) 

li) • ■ • K{ka, Is)) = (2 7 r) 3 5 2C (l 1 + . . . + la)fcz>(*i + • • • + k 2 )B 3 (h,h) (91) 

h) • ■ ■ k(&4, 1 4 )) = (27t) 3 5 2 d(1i + . . . + U)S 1D (ki + ... + k 2 )T 4 {ki, li). (92) 

The relations that we will be using most in our derivations are the orthogonality relationship of the Bessel functions Eq. dBlt and the representation of the 
2D Dirac delta function 



3i(hr\\)ji(l2r\\)dr\\ = (^^J s io(h - h); J 



e X p{irx • (li - 1 2 )} dV = (2tt) 2 6 2D (h - 1 2 ). (93) 



Our convention for the Fourier transform of an arbitrary real-space function /t(rii , rx) to Fourier space «(fc, 1) for small angular scale approximation is 
given by: 

K ( k > !) = J dr \\ J -^- k ji( kr \\) ex P(- 1 - r ±) K ( r h r -i-)' «(r|| , rx.) = J dk J ^kji{kr\\) exp(-l • r x )«(fe, 1). (94) 

We will also be working with partial transforms such as n(k, rx) and /t(rii , 1) which are defined in an obvious way. 

9.1 Flat Sky Without Mask 

In this section we will consider the power spectra associated with multi-spectra in a flat patch of sky suitable for smaller surveys. These are cross power 
spectra P pq (l), that are the Fourier transforms of cross correlation function of fields constructed from the moments of the original field «(r), e.g. K p (r) 
and n q (r). Being collapsed multipoint statistics, they carry information about the associated multispectra of order p + q, although they themselves are 
just two-point objects in real space. We will be using the same convention for the Fourier transform as the previous section: 
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K(r±) = hi «(l)exp(»lT±)d 2 l; k(1) = / k(1) exp(-il ■ r ± )d 2 r ± ; («(L>(1 2 )) C = (2tt) 2 P(0<5 2 ,d (li - 1 3 ) (95) 

In these expressions we have suppressed the explicit radial dependence - which will be introduced at the end of sections. We will first consider the skew 
spectrum. The cumulant correlator of interest in this case is (K 2 (ri)At(r2)) c . This is related to the underlying bispectrum B(l\ , l 2 , 13) and we denote the 
associated power spectrum by P21Q). 

k(2)(1)= (2^ yV(rJ]exp(-il-r ± )d 2 r ± ; K 2 (r ± ) = / K (h)K(h) exp{(li + 1 2 - 1) • r ± }d 2 lid 2 l 2 ; (96) 

k (2) (1) = ^Tf j «(li)«(la)4u>(li + I2 - l)d 2 lrd 2 l 2 . (97) 
Using the above expressions we can write down the flat-sky version of the skew spectrum: 

7>ai(I) s (k (2) (1)k*(1)) c = J B 3 {h,h,l)S(h,h,l)hhdhdh. (98) 
In deriving this we have carried out the angular integrals 4>i 1 and 0i 2 by using the following equation dHivon et alj|2002l) : 

d0i 2 &2D (li + I2 + 1) = 2tt5 (h,h,l), (99) 



where 

S(h,h, h) = -(ll + ll + ll- 2hh - 2hh - 2hh)~ 1/2 - (100) 

The derivation outlined above implicitly assumes that the bispectrum has no angular dependence in Fourier space. This is valid for the "stellar" model we 
will be considering. Carrying through the analysis in a very similar way we can write down the first of a set of two degenerate power spectra associated 
with the trispectrum P22© : 

PaaC) = <k (2) (1)« (2) (1)*)c = J (k(1iM1 2 )k(1)/<14))c<52£>(1i + h - l)S 2D (h + h + l)d\d%d 2 hd 2 h. (101) 
Using Eq.d99ll again to simplify the angular integrals we find 

V 22 (l) = J hdh J hdh J hdl 3 J hdh%(h,h,h,h)S(h,h,l)S{l 3 ,h,l). (102) 
In an analogous way, going through the same algebra for the other degenerate power spectra associated with trispectrum we find 

K (3) (l)= /"«; 3 (rpe)exp(ilT X )d 2 k; k (3) (1) = f k(1i)k(1 2 )«(13)&o(1i + h + h - l)d 2 lid 2 l 2 d 2 l 3 (103) 
P 3 i(l) = (k (3) (1)*k(1)) = J T4(h,h,h,h)S(h,h,l')S{l' ,h,l)hhhl'dhdhdhdl' . (104) 

In our derivation we have decomposed the &2D function <5 2 _d (li + 1 2 + 13 — 1) in terms of two &2D function &2D (li + 1 2 + 13 — 1) = J <5 2 r> (li + 1 2 + 
\')&2D (I3 - 1 - l')d 2 l' and used Eq.l[99j individually on each of them. We assume no angular dependence for the trispectrum, valid for the stellar model 
that we consider. For the radial dependence we need to replace the Tk(li, h, h, h) with 7~i(h, h, h, h\ rt) which will not affect the rest of the analysis. 

The power spectra Vs\(l) and V22 (I) both probe the configuration dependence of the trispectrum T4 (h ,h,h,h) in a restricted sense. While the 
power spectrum V22 (1) considers all possible configurations with the diagonal of the quadrangle (formed by momentum vectors li) constant, V31 (I) gives 
an estimate where one side of the quadrangle is kept fixed while all other sides as well as both diagonals vary. 

Given a specific model for the multispectra now we can compute the associated power spectra and compare them with simulated and observed data. 
To consider radial dependence we need to start our analysis from Eq.l|88t for the bispectrum and Eq.l|89t for the trispectrum. Following exactly similar 
analysis we will recover the flat sky versions Eq.d64t and Eq.f69l. 

9.2 Flat Sky With a Mask 

We will consider a very general mask w(r) without enforcing any symmetry. We will show that the estimated power spectra from a masked data is a 
convolved estimate of the underlying true estimates. We will derive the convolution function and how it is determined by the properties of the mask and 
develop procedures to deconvolve the effect of the mask to have an unbiased estimator. The analysis will parallel our discussion for the spherical sky. Let 
us start by defining a convolved scalar field k(r) which depends both on the original field «(r) as well as the mask w(r), «(r) = ft(r)iu(r) In the Fourier 
domain this takes the form of a convolution. If we represent the Fourier transform of k(r) by k(1) we can express it in terms of the following expression 
using the Fourier transform of the mask w(l) and the unsmoothed convergence field «(1): 

R(I)= / K(li)w(l 2 )*2D(li+l2-l)d 2 lid 2 l2, (105) 



Higher-order Statistics for 3D Weak Lensing 17 



which with the help of a kernel K\\i [w], which encodes the effect of the mask, takes a more compact form: 

«(li) = [ <fhK llh [w]; K hl2 [w]= [ d 2 hw(h)S 2D (h-h + h)- (106) 



The power spectrum associated with the masked fields is given in terms of the ker nel S(l, h, I2) and the power spectrum of the unmasked field P(h) and 
the power spectrum of the mask P W {1) — ^L J d(f>iw(l)w(l)* dHivon et alj|2002h : 

Vn(l) = ^- / #i<«(1)k*(1)) = / V(h)Pw(h)S(l,h,h)hdhhdh. (107) 



2tt 

We can rewrite the same equation in a compact form jHivon et al.ll2002h : 

fiii(h) = I V{h)M hl2 hdh; M hh =2w I P w (l)S(l,h,l 2 )ldl. (108) 



Here we have introduced the kernel Mi t i 2 and used the notation d 2 l — ldld(f> for the integration variables on the surface of the sky (2D flat patch). 
S(h, h, h) is defined in Eq. dlOOt , This result has general applicability and does not depend on any specific mask properties. Hence Pn(l) can be used 
as an estimator for the deconvolved power spectrum P(h). Carrying out a exactly similar procedure we can obtain the results for the skew spectrum 
or the kurt-spectrum. If we denote the power spectrum associated with the correlation function (k p (r\)k q (v 2 )) c by P Pq (l) and the deconvolved power 
spectrum P pq (l) which is the Fourier representation of the correlation function {K p (ri)n q (r2))c then they are related by the same expression: 

Pw('i) = ;r / d<M£(li) P £*(li) 9 } = / M hh V m [h)hdh. (109) 



2tt m 

The Fourier transforms of the squared field k 2 with a mask are related by the usual coupling matrix, 

k (2) (h) = J K hl2 [w]K i2) (h)d 2 h; « (3) (J 2 ) = J K lll2 [w]K l3) (h)d 2 h. (110) 
and masked power spectra V21H1) = (k' 2 '(1i)k(12)} = (27t) 2 <52d(1i — b) are given by 

P2i(h) = I Mi x iJ>2i{h)hdh; P 22 (h) = I M wl J>22{l2)hdh; P 3 i(h) = [ M hh V 3 i(l 2 )l2dl2. (Ill) 



These equati ons represen t the fl at-sky version of the all-sky expressions Eq.< l66t and Eq.ll70t. They generalise the results obtained for the flat-sky power 
spectrum by iHivon et al.l d2002l) . The coupling matrix Mi t i 2 introduced in this section is the flat-sky analogue of its all-sky counterpart. The power 
spectra described here, i.e. Vu(l), P21 (I), V22 (I) , V31 (I) are associated with relevant multispectra of the same order. These are useful probes of associated 
mutispectra as they do not compress the available information to a single number and retain some of the relevant shape dependence. The fact that unbiased 
estimators can be constructed by simple inversion means estimation of such multispectra from simulations and observational data may be realistically 
possible even in the presence of complicated masks with non-trivial topology. The issues of analysis of noise subtraction can be dealt with in a similar 
manner. 



10 CONCLUSIONS 

Future weak lensing surveys will play a big part in further red ucing the uncertainty in fundamental parameters, including those that describe the evolu- 
tion of equation of state of dark energy Refre gier et all feOld) . Weak lensing s urveys can exploit both the angular diameter distance and the growth o f 
structure to constrain cosmologi cal parameters, and can test the gravity model [Heavens. Kitching & Verdd ( 2007); lAmendola. Kunz & Saponej {2008); 
IBenvon. Bacn & Kovamal d2009l) . For recent results, see lSchrabback et all (2009); Kilbinge r et al.l d2009l) . Such constraints from weak lensing are com- 
plementary to those obtained from cosmic microwave background studies and from galaxy surveys as they probe structure formation in the dark sector at a 
relatively low redshift range. Initial studies in weak lensing were restricted to studying two-point functions in project ion for the entire source distribution. 
It was, however, found that binning sources in a few photometric redshift bins can imp rove the constra i ntaHul d 19991). More recently a full 3D for malism 
has been developed which uses photometric redshift of all sources without any binnin dHeavensI d2003l) ; ICastro et all d2005t) ;lHeavens et al (2006). These 
studies have demonstrated that 3 D lensing can provide more powerful and tighter constraints on the dark energy equation of state parameter, on neutrino 
masses de Bernardis et al. as well as testing braneworld and other alternative gravity models. Mos t of these 3D works have primarily focussed o n 

power spectrum analysis, but in future accurate higher-order statistic measurement should be possible (e.g. Taka da& Jainld2004l) :l Sembo loni et "ali l2009l) . 

In this paper we have generalized such studies analytically to multi-spectra which takes us beyond conventional power spectrum analysis. The 
previously obtained analytical results were developed for the statistical study of weak lensing observable using generic models for the multispectra 
of the underlying mass distribution. Later on we specialize the results for the case of specific examples using the hierarchical ansatz, where higher- 
order multispectra are constructed from various products of power spectra organized in all possible topological diagrams with different amplitudes. The 
analytical results are developed both for near all-sky surveys as well as for flat patches of the sky. The formalism developed does not depend on the 
background cosmology and can be used to predict level of non-Gaussianity for both primary as well as secondary non-Gaussianity. 

The higher-order multispectra contain a wealth of information in through their shape dependence. Though partly degenerate, this information can be 
invaluable for constraining structure formation scenarios. However determination of the multispectra and their complete shape dependence is not an easy 
task from noisy data. In this paper we advocate a set of statistics called "cumulant correlators" which were first used in real space in the context of galaxy 
surveys and later extended to CMB studies. Here we have presented a general formalism for the study of the power spectra or the Fourier transforms of 
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these correlators. We present a 3D analysis which takes into account the radial as well as on the surface of the sky decomposition. We start by relating 
various representation of multispectra in three dimensions. We relate the spherical representation and the Fourier representations with other possibilities: 
mixed modes of representations. These allow us to relate the harmonic decomposition of convergence directly with that of underlying mass distribution. 

We have restricted this study to the third and fourth order, though it can be generalised to higher order and some of our results are valid at arbitrary 
order. At third order, we define a power spectrum which compresses information associated with a bispectrum to a power spectrum. This power spec- 
trum Cf' 1 (r2, ri) is the cross-power spectrum associated with squared convergence maps K 2 (ri, Cl) constructed at a specific radial distance n against 
K 2 (r2, Cl) at T2 In a similar manner we also associate power spectra C 2 ' 2 and Cf' 1 with associated trispectra T/ 1 /* (L; n). There are two different power 

spectra at the level of trispectra which are related to the respective real-space correlation functions (n 2 (ri, Cl)n 2 (r2, Cl')} and (ft 3 (ri, 0)fv(r2, "))• 
We expressed these real-space correlators in terms of their Fourier space analogue which take the form of C, 3 ' 1 ^, ri) and C 2 ' 2 (r2, ri). We develop 
analytical expressions to take into account the photometric redshift errors in these power spectra. While we present formalisms which are completely 
general, we also use the Limber approximation to reduce the dimensionality of the relevant integrations. These when combined with specific hierarchical 
models of gravitational clustering can make analytical results remarkably simpler. 

The statistics presented here will be a useful tool in studying non-Gaussianity in alternative theories of gravity; which are one of the important 
science drivers for the future generations of weak lensing surveys. We plan to present detailed results elsewhere in future. 
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APPENDIX A: REALISTIC SELECTION FUNCTION AND PHOTOMETRIC REDSHIFT ERRORS 

The results obtained in the main text was simplified for clarity, ignoring the fact that in a realistic survey, the average number density of sources will 
decline with distance, and the distances estimated from photometry will include errors. We consider these here. 

Al All Sky results 

The lensing potential can only be sampled at the position of galaxies. Hence it can be written as sum over galaxy positions. This discrete sum can be 
expressed as 



Here W is an arbitrary weight function, and r g is the distance to galaxy g assuming a fiducial cosmology. The convergence depends of course on teh 
correct distance in the true cosmology, r. Replacing the discrete sum over the galaxy positions with an integral we can write 




(Al) 



9 




(A2) 
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The quantity n(r) comprised of sum of delta functions which peaks at observed positions of the galaxies. Ensemble averaging of these quantities will 
reduce the equation to 

(A3) 



«im(*;r) = y- / d A r"n{r)K{r)k3i{kr u )Y lm (Sl g )W{r") 

Because of the discrete nature of source galaxies the estimator will have a scatter due to the shot noise. It will also have contribution from source 
clustering. While we will investigate the effects of photometric redshift errors, we will ignore the uncertainties in the photometric redshift distribution of 
sources which means we can write n(r°)d 3 r° = n z (z p )dz p dfl / 'A-k . We will ignore the effect of source clustering which does not play a dominant role 
in the error budget for deep surveys. 

Kp m (fc;r) = j dz d Q n z (z p )K(r)kji(kr°)Yi m (£l)W(z p ) (A4) 

where r° is the fiducial distance at redshift z p . 

The primary effect of photometric redshifts is to smooth the source distributions along the line of sight distribution. If p(z\z p ) denotes the probability 
of the true redshift being z given the photometric redshift z p , the above equation, when modified to take into account the effect of photometric redshift 
error can be written as: 



K? m (k;r) — y' J dz p J dz J dttn z (z p )p(z\z p )>t(r)kji(kr )Yi m (&)W(z p ) 

After expanding the ft(r) and carrying out the angular integrations we can eventually arrive at the following expression: 

K? m {k;r) = J / dz p / dzn z (z p )p(z\z p )kji(kr°) / dk'k'ji(k'r°)W(z p )Ki m (k';r) 



(A5) 



(A6) 



Typically p(z\z p ) is modelled as a Gaussian for simplicity, though it may have catastrophic failures. These can be included by modification of p(z\z p ) 



p(z\z p ) = 



2-ku z (z 



exp 



(z p Z -(- Zbi as ^ 

2a 2 z (z) 



(A7) 



In this expression zuas is the possible bias in the photometric redshift calibration. The dispersion in error a z (z) depends on the redshift. The photometric 
redshift errors evidently introduces error in the radial direction. Having expressed the nf m (k; r) by taking into account the photometric redshift errors in 
terms of K; m (fe, r), we can now construct the multispectra for the observed harmonics by relating Ki m (k; r) to Si m (k; r) as outlined in the main text . 



A2 Flat Sky Expressions 

Here we give fiat-sky results. We start by decomposing the convergence field n(k, 1) as before. Here 1) is the observed convergence assuming a 

fiducial cosmology. 



«°(Jfe,l) = 



2tt 



dz p n z (z p )n(v)kji{kr°) exp[— il ■ 9]. 



By expressing k{t\\ , rj_) in terms of n(k, 1) we find 
dn f . f dk' 



K°(k,l) = - 

7T 



K°(k,l)=- 
7T 



dz,- 



kji {kr || )k'ji (fc'r || )n z (z p )K,(k' , 1) 



dz / dz p p(z p \z)n z {z p )k ji(kr\\°) 



dm 



" k -.k' ji(k'r')K(k' , 1). 



2tt 



(A8) 



(A9) 



(A10) 



This expression enables us to relate the theoretical predictions for a fiducial cosmology to the observed k° (k, 1). We have assumed n(r)drj_d 2 r^ — 
\n z (z p )dz p d? V||, where A is the solid angle of the sky covered. The statistical properties such as the bispectrum and trispectrum of the field n(k, 1) 
derived in the main text can now be used to predict the observed statistics of no(k, 1). Mixing of modes due to the photometric redshift error will couple 
the radial modes, whereas the partial sky coverage mixes angular modes. 



APPENDIX B: USEFUL MATHEMATICAL RELATIONS 
Bl Spherical Bessel Functions 

The orthogonality relationship for the spherical Bessel functions is given by the following expression: 



k 2 ji(kn)ji(kr2)dk = 



5id(ti — r 2 ). 



2(1 + 1/2) 2 _ 

The extended Limber approximation is also implemented through the following approximate relation lLo Verde & Afshordil (2008): 



(Bl) 



/ 



F(k)ji(kri)ji(kr 2 )dk 



2rj 



F (^ — ^j 5 1D (n - r 2 ). 
Thus for high I the spherical Bessel functions can be replaced by a Dirac delta function 8m ■ 

}TJ l{x) = \f^i 6lD { l + \- x )- 

B2 Spherical Harmonics 

The completeness relationship for the spherical harmonics is given by: 

Y lm (Cl)Y lm {Cl') = 8 2D ((l - Cl'). 

lm 

The orthogonality relationship is as follows: 
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(B2) 



/ 



dClY lrn (Cl)Y l , m ,(Q) = 8^8, i 



K £ K 
mm' 



B3 3J Symbols 

The following properties of 3 J symbols were used to simplify various expressions. 



V (2/3 + 1) ( h h h )( h , h , 1 ) 
t-^i y mi m 2 vnz J y m 1 m 2 m J 



trK trK 



h m 3 

E 

m\m2 

(-1)" 



h h I3 \ I h h I3 \ h l ' 3 m : 



mi m 2 m:i 

I I I' 

rn —m 



mi m 2 m 3 



(-1)' ,K 



2l A + l 



(B3) 



(B4) 



(B5) 



(B6) 
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